/*
 * EquilibriumSolver2.cpp
 *
 *  Created on: 26 Jul 2011
 *      Author: Allan
 */

#include "EquilibriumSolver2.h"

// C++ includes
#include <cmath>
#include <functional>
#include <iomanip>
#include <memory>
using namespace std;
using namespace std::placeholders;

// Boost includes
#include <boost/foreach.hpp>
using namespace boost;

// GeoReact includes
#include "Utils/Algorithms.h"
#include "Utils/SpeciesMap.h"
#include "Utils/SpeciesUtils.h"

EquilibriumSolver2::EquilibriumSolver2() :
tolerance(1.0E-6), maxiterations(100)
{}

EquilibriumSolver2::EquilibriumSolver2(const Multiphase& multiphase, const ReactionSystem& reactions) :
tolerance(1.0E-6), maxiterations(100)
{
	Initialiaze(multiphase, reactions);
}

void EquilibriumSolver2::Initialiaze(const Multiphase& multiphase, const ReactionSystem& eReactions)
{
	// Set data members
	this->multiphase = multiphase;
	this->eReactions = eReactions;
	
	// Initialiaze the stoichiometric matrix of the system of equilibrium-controlled reactions
	uei = eReactions.AssembleStoichiometricMatrix(multiphase.GetSpecies());
	
	// Initialiaze the Newton solver
	newtonsolver.Initialiaze(multiphase.GetNumSpecies());
}

void EquilibriumSolver2::SetDecompositionAlgorithm(NewtonSolver::Decomposition decomposition)
{
	newtonsolver.SetDecompositionAlgorithm(decomposition);
}
	
void EquilibriumSolver2::SetConvergenceOptions(double tolerance, uint maxiterations)
{
	this->tolerance = tolerance;
	this->maxiterations = maxiterations;
}

void AppplyUnderRelaxation(const VectorXd& nOld, VectorXd& n)
{
	for(unsigned i = 0; i < unsigned(n.rows()); ++i)
	{
		const double delta  = n[i] - nOld[i];
		const double factor = max(1.0, -2.0 * delta/nOld[i]);
		
		n[i] = nOld[i] + delta/factor;
	}
}

ConvergenceDiagnostics EquilibriumSolver2::Solve(const EquilibriumConditions& econditions, double T, double P, VectorXd& n)
{
	// Auxiliary variables 
	VectorXd nOld, a, ln_a, ln_Ke, hj;
	
	// Calculate the natural logarithm of the equilibrium constants of the equilibrium-controlled reactions
	ln_Ke = eReactions.EquilibriumConstants(T, P).array().log();
	
	// The equilibrium condition matrix
	const MatrixXd& alpha = econditions.GetConditionMatrix();
	
	// The current state of the convergence
	ConvergenceDiagnostics diagnostics = {0, false};
	
	while(!diagnostics.converged && diagnostics.iterations++ < maxiterations)
	{
		// Calculate the activities of all species
		a = multiphase.Activities(T, P, n);
		
		// Calculate the natural logarithm of the activities
		ln_a = a.array().log();
		
		// Calculate the (hj) vector with the equilibrium condition functions
		hj = econditions.EvaluateConditionFunctions(T, P, n, a);
		
		// Create an equilibrium problem instance to calculate the function F and jacobian J
		EquilibriumProblem2 problem(uei, alpha, hj, ln_a, ln_Ke);
		
		// Update the previous composition state (nOld)
		nOld = n;
		
		// Perform one iteration with the Newton's method and retrieve the convergence state
		diagnostics.converged = newtonsolver.Iterate(problem, n);
		
		// Ensure non-negative values of compositions
		AppplyUnderRelaxation(nOld, n);
	}
	
	return diagnostics;
}

